% read files and create ETData Objects
if ~exist('e24', 'var')
    % change workfolder
    cd('D:\\Workspace\\Compensation\\Data\\CK50S');

    % create object
    e24 = ETData('ST20181124W');
    e26 = ETData('ST20181126WC');
    e27 = ETData('ST20181127WC');
    e28 = ETData('ST20181128W');

    % add compensation
    e26.comp([1 2 3]); e27.comp([1 2 3]);

    f1 = e24.plot; f2 = e26.plot; f3 = e27.plot; f4 = e28.plot;
end

% model
if ~exist('B', 'var') || isempty(B)
    sx = [1 2 4 6 7];
    % sx = [3 4 6 7 8];
    sy = 1;

    [X1, Y1, U1, V1] = model(e24, [1, 22]);
    [X2, Y2, U2, V2] = model(e26, [1, 17]);
    [X3, Y3, U3, V3] = model(e27, [1, 17]);
    [X4, Y4, U4, V4] = model(e28, [1, 14, 25]);
    % X1 = e24.Temp(:, sx); Y1 = e24.Eddy(:, 1);
    % X2 = e26.Temp(:, sx); Y2 = e26.Eddy(:, 1);
    % X3 = e27.Temp(:, sx); Y3 = e27.Eddy(:, 1);
    % X4 = e28.Temp(:, sx); Y4 = e28.Eddy(:, 1);

    [X, Y   ] = modelxy(X1, Y1, X2, Y2);
    [X, Y, B] = modelxy(X, Y, X3, Y3);
    [X, Y, B] = modelxy(X, Y, X4, Y4);

    U = [U1; U2; U3; U4];
    V = [V1; V2; V3; V4];
    D = regress(V, U);

    B(1) = [];
end

% plot
if ~exist('l1', 'var')
    figure(f1); subplot(2, 1, 1);
    x  = e24.Temp(:, sx);
    c  = e24.Eddy(1, 1) - x(1, :) * B;
    y  = x * B + c;
    dt = (x(22, :) - x(21, :)) * (D - B);
    y(22:end) = y(22:end) + dt;
    disp(max(abs(e24.Eddy(:, 1) - y)));
    l1 = plot(e24.Time, y, '--');

    figure(f2); subplot(2, 1, 1);
    x  = e26.Temp(:, sx);
    c  = e26.Eddy(1, 1) - x(1, :) * B;
    y  = x * B + c;
    dt = (x(17, :) - x(16, :)) * (D - B);
    y(17:end) = y(17:end) + dt;
    disp(max(abs(e26.Eddy(:, 1) - y)));
    l2 = plot(e26.Time, y, '--');

    figure(f3); subplot(2, 1, 1);
    x  = e27.Temp(:, sx);
    c  = e27.Eddy(1, 1) - x(1, :) * B;
    y  = x * B + c;
    dt = (x(17, :) - x(16, :)) * (D - B);
    y(17:end) = y(17:end) + dt;
    disp(max(abs(e27.Eddy(:, 1) - y)));
    l3 = plot(e27.Time, y, '--');

    figure(f4); subplot(2, 1, 1);
    x  = e28.Temp(:, sx);
    c  = e28.Eddy(1, 1) - x(1, :) * B;
    y  = x * B + c;
    dt = (x(14, :) - x(13, :)) * (D - B);
    y(14:end) = y(14:end) + dt;
    dt = (x(25, :) - x(24, :)) * (D - B);
    y(25:end) = y(25:end) + dt;
    l4 = plot(e28.Time, y, '--');
end
